#work out whether SNPs dn or ds
#reworked version to use Biostrings not GeneR

#wrapper for a per contig analysis
codingSNPref <- function(AllSNPref,dir.embl="/Paterson/Datafiles/grouse/embl_exon2",seq.name=NULL,file.name=NULL){
	#AllSNPref is an object (dataframe) from baseCoverage
	#expects embl files in dir.embl, perhaps created by add.CDS
	
	if(is.null(AllSNPref)) return('null.input')
	if(identical(AllSNPref,0)) return('zero.input')
	
	
	if(is.null(seq.name)) seq.name <- AllSNPref[1,"cDNA"]
	if(is.null(file.name)) {
		#fasta.file.name<-paste(dir.fasta,"/",seq.name,"_matched.fasta",sep="")
		embl.file.name<-paste(dir.embl,"/",seq.name,"_exon.embl",sep="")
		}else{
		#fasta.file.name <- file.name[1]
		embl.file.name <- file.name[1]
		}
	
	if(!file.exists(embl.file.name)){
		warning(embl.file.name, " does not exist... returning\n")
		AllSNPref$cds <- FALSE
		AllSNPref$dn <- FALSE 
		return(AllSNPref)
		}

	#require(Biostrings)
	
	#for some reason won't read straight into buffer
	
	#ref.seq <- readEmblSeq(embl.file.name)
	#ref.seq.set <- read.DNAStringSet(fasta.file.name)
	#ref.seq <- DNAString(ref.seq)
	
	#extract feature table
	tmp.features <- system(paste("grep \'FT\'",embl.file.name,"|sed \'s/FT[ \t]//\'"),intern=TRUE)
	cds.lines <- grep('^[ ]*CDS ',tmp.features)
	product.lines <- grep('/product',tmp.features)
	grep.cds <- character(length=length(cds.lines))
	for(cdsi in 1:length(grep.cds)){
		grep.cds[cdsi] <- paste(tmp.features[cds.lines[cdsi]:(product.lines[cdsi]-1)],collapse="")
		
		}
	#extract coding seqs
	grep.cds <- sub('^[ ]*CDS ','',grep.cds)
	grep.cds <- gsub('[ ]*','',grep.cds)

	#grep.cds <- system(paste("grep \'CDS\'",embl.file.name,"|sed \'s/FT[ \t]*CDS[ \t]*//\'"),intern=TRUE)
	complementCDS <- FALSE
	if(substr(grep.cds[1],1,4)=="comp"){
		complementCDS <- TRUE
		}

	ref.seq <- readEmblSeq(embl.file.name)
	
	grep.cds <- sub("complement\\(","",grep.cds)
	grep.cds <- sub("join\\(","",grep.cds)
	grep.cds <- gsub(")","",grep.cds)

	#find masks
	mask.lines <- grep("/note=\"mask\"",tmp.features) -1
	if(length(mask.lines)>0){
		grep.mask <- tmp.features[mask.lines]
		grep.mask <- sub('^[ ]*misc_feature','',grep.mask)
		grep.mask <- gsub('[ ]*','',grep.mask)
		
		mask.list <- strsplit(grep.mask,"\\.\\.")
		mask.list <- lapply(mask.list,as.numeric)
		mask.list <- lapply(X=mask.list,FUN=function(X){
			if(length(X)==1){
				return(c(X,X))
				}else{
				return(X)
				}
			})
		mask.list <- matrix(unlist(mask.list),ncol=2,byrow=T)
		colnames(mask.list) <- c("start","end")
		mask.list <- as.data.frame(mask.list) 
		#output starts and ends of mask as a dataframe
		tmp.loc <- numeric(0)
		for(i in 1:nrow(mask.list)) tmp.loc <- c(tmp.loc,seq(mask.list[i,1],mask.list[i,2]))
		for(i in 1:length(tmp.loc)){
			ref.seq <- paste(substr(ref.seq,1,(tmp.loc[i]-1)),"#",substr(ref.seq,(tmp.loc[i]+1),nchar(ref.seq)),sep="")
			#will break if masked nucl at end of sequence
			}

		}



	#replace masked nucleotides with a "#"
	#if(length(tmp.loc)>0){
	#			}

#	AllSNPref$cds <- FALSE
#	AllSNPref$dn <- FALSE
	#AllSNPref[,c('cds','dn')] <- t(sapply(X=AllSNPref[,"id"],FUN= mutateSeq,cds.str=grep.cds,ref.seq=ref.seq,complementCDS=complementCDS))
	tmp <- t(sapply(X=AllSNPref[,"id"],FUN= mutateSeq,cds.str=grep.cds,ref.seq=ref.seq,complementCDS=complementCDS))
	tmp <- as.data.frame(tmp)
	names(tmp) <- c("cds","dn")
	cbind(AllSNPref,tmp)
	}



#seq.name <- "cDNA_260-1"
#dir.embl <- "/Paterson/Datafiles/grouse/embl_exon"

#read sequence into buffer




#sequence 0 for original seq, seq 1 for mutated seq, seq 2 for temporary seq
#translateCDS <- function(cds.str,seqno=0){
	#expects cds.str in for 1..10,20..30 etc
	#and for a sequence in buffer seqno
#	tmp.from <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('\\.\\.[0-9]*$','',X)}))
#	tmp.to <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('^[0-9]*\\.\\.','',X)}))
#	assemble(seqno=seqno,from=tmp.from,to=tmp.to,strand=0,destSeqno=2)
#	translate(2)
#	}
	
mutateSeq <- function(snp.id,cds.str,ref.seq,complementCDS=FALSE){
	#expects a snp.id, ie id from baseCoverage
	#cat(snp.id,"  ")
	require(Biostrings)
	snp.id <- sub("^cDNA_","cDNA--",snp.id)
	snp.id <- sub("^ssh_","cDNA--",snp.id)
	tmp.from <- list()
	tmp.to <- list()
	for(i in 1:length(cds.str)){
		tmp.from[[i]] <- as.numeric(sapply(X=strsplit(cds.str[i],","),FUN=function(X){sub('\\.\\.[0-9]*$','',X)}))
		tmp.to[[i]] <- as.numeric(sapply(X=strsplit(cds.str[i],","),FUN=function(X){sub('^[0-9]*\\.\\.','',X)}))	
		}
	
	
	tmp.id <- strsplit(snp.id,"_")[[1]]
	#if(substr(tmp.id[1],2,6)=="ontig"){
	#	mut.start <- as.numeric(tmp.id[2])
	#	}else{
	mut.start <- as.numeric(tmp.id[2])
	#	}
	
	
	mut.cds <- FALSE #test whether in coding region
	mut.dn <- FALSE
	for(i in 1:length(tmp.from)){
		if(any(mut.start>tmp.from[[i]]&mut.start<tmp.to[[i]])) mut.cds <- TRUE
		}

	if(mut.cds==FALSE) return(c(as.logical(mut.cds),as.logical(mut.dn)))
	#if the mutation isn't in a coding region just return
	
	
	#mutation in coding region...
	mut.old <- sub("-","",tmp.id[5])
	mut.new <- sub("-","",tmp.id[6])
	if(nchar(mut.old)!=nchar(mut.new)){
		#if there's an indel in the coding sequence assume non-synonymous and return
		mut.dn <- TRUE
		return(c(as.logical(mut.cds),as.logical(mut.dn)))
		}

	#should now be carrying on only if a straight swap of nucleotide(s)
	mut.len <- nchar(mut.old)
	mut.stop <- mut.start + mut.len -1

	#will probably fail for large sequences
	cds.old <- character(length=length(tmp.to))
	
	#must be quicker way, but generally only 1 to 3 exons
	#translate ref sequence
	cds.old <- character(length=length(tmp.to))
	for(cds in 1:length(tmp.to)){
		for(j in 1:length(tmp.to[[cds]])){
			cds.old[cds] <- paste(cds.old[cds],substr(ref.seq,start=tmp.from[[cds]][j],stop=tmp.to[[cds]][j]),sep="")
			}
		} 
		

	cds.old <- gsub("#","",cds.old)
	
	if(complementCDS){
		cds.old <- as.character(reverseComplement(DNAStringSet(cds.old)))
		}
	
	
	#translate, but use GeneR function to allow Ns
	old.prot <- as.character(sapply(cds.old,GeneR::strTranslate))
	
	new.seq <- paste(substr(ref.seq,1,(mut.start-1)),mut.new,substr(ref.seq,(mut.start+mut.len),nchar(ref.seq)),sep="")

	#new.seq <- replaceLetterAt(ref.seq,at=seq(mut.start,mut.start+mut.len-1),letter=mut.new)
	
	#translate mutated sequence
	cds.new <- character(length=length(tmp.to))
	for(cds in 1:length(tmp.to)){
		for(j in 1:length(tmp.to[[cds]])){
			cds.new[cds] <- paste(cds.new[cds],substr(new.seq,start=tmp.from[[cds]][j],stop=tmp.to[[cds]][j]),sep="")
			}
		} 
		

	cds.new <- gsub("#","",cds.new)
	
	if(complementCDS){
		cds.new <- as.character(reverseComplement(DNAStringSet(cds.new)))
		}
	
	
	#translate, but use GeneR function to allow Ns
	new.prot <- as.character(sapply(cds.new,GeneR::strTranslate))
	
	if(!identical(old.prot,new.prot)) mut.dn <- TRUE
	
	c(as.logical(mut.cds),as.logical(mut.dn))
	}


catagoriseSNPs <- function(snp.id,grep.cds,ref.seq,complementCDS=FALSE){
	tmp.mut.list <- vector(mode="list",length=length(grep.cds))
	for(i in 1:length(grep.cds)){
		tmp.mut.list[[i]] <- mutateSeq(snp.id,grep.cds[i],ref.seq,complementCDS=complementCDS)
		}
	c(any(sapply(tmp.mut.list,function(X){X$mut.cds})),
		any(sapply(tmp.mut.list,function(X){X$mut.dn})))

	}


readEmblSeq <- function(file.name){
	tmp.lines <- readLines(file.name)
	tmp.lines2 <- tmp.lines[(grep('^SQ',tmp.lines)+1):(grep('//',tmp.lines)-1)]
	rm(tmp.lines)
	tmp.lines2 <- gsub(' ','',tmp.lines2)
	tmp.lines2 <- gsub('[0-9]','',tmp.lines2)
	paste(tmp.lines2,collapse="")
	}


